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Wc use a separable mode expansion estimator with WMAP data to estimate the bispectrum for aU 
the primary famihes of non-Gaussian models, including non-scaling feature models, the flattened 
model, DBI and ghost inflation, as well as previously constrained simple cases. We review the late- 
time mode expansion estimator methodology which can be applied to any non-separable primordial 
and CMB bispectrum model, and we demonstrate how the method can be used to reconstruct the 
CMB bispectrum from an observational map. Wc extend the previous validation of the general 
estimator using local map simulations. We apply the estimator to the coadded WMAP 5-year V 
and W channel maps, reconstructing the WMAP bispectrum using I < 500 multipoles and n — 31 
orthonormal 3D eigenmodes; both the mode expansion parameters and the reconstructed 3D WMAP 
bispectrum are plotted. We constrain all popular nearly scale-invariant models, ensuring that the 
theoretical bispectrum is well-described by a convergent mode expansion. Constraints from the local 
model /nl = 54.4 ± 29.4 and the equilateral model /nl = 143.5 ± 151.2 ( Fnl = 25.1 ± 26.4) are 
consistent with previously published results. (Here, we use a nonlinearity parameter i%L normalised 
to the local case, to allow more direct comparison between different models.) Notable new constraints 
from our method include those for the constant model Fnl = 35.1 ± 27.4, the flattened model 
^NL = 35.4 ± 29.2, and warm inflation Fnl = 10.3 ± 27.2. We investigate feature models, which 
break scale invariance, surveying a wide parameter range for both the scale and phase, and we 
flnd no significant evidence of non-Gaussianity for all cases well-described by the given eigenmodes 
(i.e. with features 1* > 150). Wc propose a measure Fnl for the total integrated bispectrum and 
find that the measured value is consistent with the null hypothesis that CMB anisotropics obey 
Gaussian statistics. We argue that this general bispectrum survey with the WMAP data represents 
the best evidence for Gaussianity to date and we discuss future prospects with higher precision and 
resolution, notably from the Planck satellite. 

I. INTRODUCTION 

In an earlier paper [1] , we described a general approach to the estimation of non-separable CMB bispectra 

using separable mode expansions. Om aim here is to directly estimate the full CMB bispectrum from WMAP 
data, to survey and constrain current non-Gaussian primordial theories, and to discuss the prospects for 
reconstructing the bispectrum with forthcoming data, such as the Planck experiment. Previous work by 
other groups has endeavoured to measure the bispectrum by using specific estimators tailored to investigate 
particular separable models, such as the well-known local and equilateral bispectra. This restriction to 
separable cases was for reasons of calculational simplicity to make the data analysis tractable, that is, 
reducing it from 0{l^^^ to 0(/max) operations. We summarise constraints that have been obtained to date 
using these methods later in section V, when we survey theoretical models; it is sufficient at this point 
to note that the present WMAP7 constraint — 10 < /nl < 74 [2] (95% confidence) does not provide any 
significant evidence for a primordial local bispectrum signal, and nor do constraints on the few other models 
investigated to date (sec the review [3]). 

Two significant developments mean that we can move beyond these specific estimators and consider 
a more general approach which includes the reconstruction of the whole bispectrum directly from the 
observational data. First, explicit calculations of the reduced CMB bispectrum hi^i^i^ in a wide-ranging 
survey of primordial theories [4, 5], demonstrated that the resulting coherent patterns of acoustic peaks 
could be represented by rapidly convergent mode expansions with a limited number of terms (irrespective of 
whether the primordial bispectrum was separable). Secondly, these complete orthonormal mode expansions 
could be transformed into a non-orthogonal frame with separable basis functions [1] in which the same 
simplifications could be exploited to efficiently calculate the estimator (11) in ©(/fnax) operations, again for 
arbitrary non-separable theoretical bispectra bij^i^. We shall employ this mode expansion methodology in 
this paper, convolving observational maps with the separable basis functions and then reconstructing the 
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observed bispectrum h^i^i,^ in an expansion using the resulting mode coefficients. Rather than looking in 
just a few specific directions within the large space of possible bispectra, this general mode decomposition 
encompasses all bispectra up to a given resolution. Our aim is to determine whether there is evidence for 
any bispectrum causing a departure from Gaussianity in the WMAP data. Of course, wc can compare with 
previous constraints for the local and equilateral models, but an important byproduct is a set of entirely 
new constraints on a wide range of non-separable models. 

While we believe this work represents a significant step forward, we note that this analysis is far from the 
last word on CMB non- Gaussianity, not least because much higher quality and higher resolution data will 
soon be available from Planck. We also note that we have only used WMAP5 data out to 1=500, together 
with a pseudo-optimal analysis of the noise and masking contributions. This paper should be considered 
primarily as a proof of concept implementation of these methods, leaving up to an order of magnitude 
discovery potential available for bispectrum signals with new CMB data, let alone future galaxy and other 
3D surveys where this approach can also be applied. We note that there are other recent methodologies in 
the literature which, in principle, can be used to extract information from the bispectrum beyond simple 
separable cases, including the bispectrum power approach of ref. [6]), bispectrum binning used in ref. [7] 
and wavelet approaches (see the review [3]). 

In section II we review general results regarding primordial and angular bispectra and their optimal esti- 
mation. The eigenmodc decomposition of the bispectrum that constitutes the foundation of our methodology 
is summarized in section III. Wc then show in section IV how this expansion can be used to reconstruct the 
full bispectrum from the data, and proceed to validate this result in section V, before directly extracting 
the bispectrum from WMAP data in section VI. We finally turn our attention to estimates of /nl for a 
wide variety of shapes, including both scale invariant bispectra (section VII) and scale-dependent oscillatory 
bispectra (section VIII). Finally, before drawing our conclusions in section X, we discuss in section IX a 
possible way to use our mode expansion technique to define a model independent constraint on the total 
integrated bispectrum extracted from the data. 



II. CMB BISPECTRUM ESTIMATION 
A. Primordial and CMB bispectrum 

TempcratTire anisotropies are represented using the aim coefficients of a spherical harmonic decomposition 
of the cosmic microwave sky, 

AT 

-jri^) = Z^aimYim{n) , (1) 

Im 

with an (ideal) angular power spectrum C; = Ylm^l'mO'l-m- The CMB bispectrum is the three point 
correlator of the a;^, 

^mirn2m3 ~ (^hm-iO'hm^'^hmi > (2) 

where, here, we assume that the -B^'j^^mama coefficients are not an ensemble average but, instead, directly 
calculated using the a/m's from a high resolution map (or maps), that is, from an experiment such as 
WMAP or Planck. We shall assume for the moment that if there is a non-trivial bispectrum then it has 
arisen through a physical process which is statistically isotropic, so we can employ the angle-averaged 
bispectrum Bij^i^ without loss of information, that is [8], 
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where hi^^i^i^ is a geometrical factor, 



Khh = \l ^ ( ) ' ^^■> 



and Qrii\mX-, is ^he Gaunt integral, 



'mim2ms 



= J dn (n) y^^^^ (n) Yi^msir^) 



- Khh ( ^\ j > (5) 

with the usual Wigncr-3j symbol (mimsma). It is more convenient to eliminate the geometrical factors 
entirely and to work with the reduced bispectrum which is defined as 

Khh = ^Z^is-^'i's's • (6) 

It is important to note the relationship between the late-time CMB bispectrum bij^i,^ and the primordial 
bispectrum B^(ki, k2, k^) from which it would arise in many models, notably inflation. The convention has 
been to remove a k~^ scaling by defining a shape function: 

S{ki,k2,k3) = ^{kik2k3fB^{ki,k2,k3). (7) 

The shape function (7) is particularly pertinent for scale-invariant models because their momentum depen- 
dence is restricted entirely to planes transverse to the diagonal k = ^{ki -1-^2-1- ^3). The CMB bispectrum 
induced by the primordial shape S is obtained from the convolution [9]: 

kihh = (^^^ J x^dx J dkidk2dk3S{ki,k2,k3) 

X A;^ (ki) Ai^ {k2) (/cs) ji^ {kix) ji^ {k2x) ji, {k^x) , (8) 

where A/(fe) is the transfer function. 

The impact of the transfer functions in (8) is to impose a series of acoustic peaks on the underlying 
primordial shape, as illustrated for the CMB bispectrum of the constant model S{ki,k2,k3) = 1 in fig. 1. 
Here, we can observe a large primary peak when all the li ~ 220. In principle, the CMB bispectrum 
is difficult to evaluate since (8) represents a four-dimensional integral over highly oscillatory functions. 
However, the integral breaks down into a product of one-dimensional integrals if the shape function is 
separable, that is, if it can be represented in the form S{ki,k2,k3) = X {ki)Y {k2) Z (k^) . In the large-angle 
limit with Ai(k) = {I <C 200) it is possible in some separable models to obtain analytic solutions, such 

as that for the constant model [5] 

A2 
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^/i/2«3 277V {2h + 1)(2£2 + 1)(24 + 1) [h+h+h + ^^ h\ ' 

This particular regular solution is important because we divide by it when plotting the CMB bispectrum 
^lihh/^hkh^^"^^ throughout this paper. Normalising with the constant model (9) is analogous to multiplying 
the power spectrum Ci's by l{l + 1), because it serves to remove an overall fr^ scaling for all scale-invariant 
bispectra, preserving the effects of the oscillating transfer functions without introducing spurious transverse 
momentum dependence. 



B. CMB bispectrum estimators 



Now it is usually presumed that the full bispectrum for a high resolution map cannot be evaluated 
explicitly because of the sheer number of operations involved O(^max)' well as the fact that the signal 
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Figure 1: The reduced CMB bispectrum for the constant model b'[°l^^l^ arising from the convolution of the primordial shape 

function S{ki,k2,ks) = 1 with transfer functions (normalised relative to the large-angle constant solution ^^"("jg''''' given in 
(9)). On the left, the 3D bispectrum is plotted over the allowed tetrahedral region of multipole triples (see fig. 2) using several 
density contours (light blue positive and magenta negative) out to h < 2000. On the right, a transverse triangular slices 
through the bispectrum is shown for h + h + h = 4000 (Planck resolution). Note the coherent pattern of acoustic peaks with 
a dominant primary peak in a broad diagonal region around h = h = h = 220. This constant model bispectrum plotted is the 
analogue of the angular power spectrum Ci 's for a purely scale-invariant model. 



will be too weak to measure individual multipoles with any significance. Instead, we essentially use a least 
squares fit to compare the bispectrum of the observed aim^s (1) with a particular (separable) theoretical 
bispectrum bf^i^^^, 



Here, bf^i^i^ will be recovered as the expectation value from an ensemble average over realisations 
or simulations created with the given reduced bispectrum. Formally, taking into account the fact that 
instrument noise and masking can break rotational invariance, the result is the general optimal estimator 
[10-12] 



li,ini 



where is the inverse of the covariance matrix Ci-^rnx,i2m2 = (f^^^imi 0/2^2) ^-iid is a suitable normal- 
isation (discussed further below). Here, we follow ref. [13, 14], by assuming a nearly diagonal covariance 
matrix {Ci^rni,i2m2 ^ Ci5i^i^ 6mi-m2) aiid approximating the estimator (11) as 

~ C C C \"'hmiO'l2m20'hm3 ~ O O;^^^ 0/3^3 ; , (1^) 
liirii '1 '2 '3 

where the tilde denotes the modification of C; and h^i^i^ to incorporate instrument beam and noise effects 
through 

Ci = bfCi + Ni and kj^i^ = hfy^hi^ k^i^i^ . (13) 
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For a relatively small galactic mask (leaving a large fraction /gky of the full sky), it has also been shown to 
be a good approximation to renormalisc using 

Khk = fs^^yblM. and = /,kyQ • (14) 

(We shall assume noise, beam and mask inclusion henceforth and drop any special notation.) Here, the 
second linear term in (12) ensures subtraction of spurious inhomogeneous noise and masking contributions 
by using the covariance matrix Cf™^ ^^^^ from an ensemble average of Gaussian maps in which these effects 
are incorporated. 

If the theoretical bispectrum bf^i^i^ has the property of primordial separability then it has been noted 
that the summation in (12) becomes much more tractable taking only ©(/^ax) operations [10]. Essentially 
this exploits the separability of the Gaunt integral (5), as well as primordial counterparts, to reduce the 
dimensionality of the integrals and summations involved in evaluating (12) (see ref. [1] for a more detailed 
discussion on this point). To date, such separability has been a property of all the primordial theories 
constrained observationally with most attention given to the canonical local model. 



Jnl normalisation 



It remains to briefly discuss the normalisation factor in (11). In the past this has been taken on a 
case-by-case manner for a given theoretical bispectrum h/J^^^ia *° 



/j2 ^th 2 

As we discuss below, this has yielded very model-dependent results for the measurement of the nonlinearity 
parameter = £. Instead, we have proposed the parameter J^nl which is much easier to compare between 
models, because it measures the integrated CMB bispectrum signal relative to that from the canonical local 
model with /nl = 1- In this case, we define[l] 

F^^ = £, with Ar2 = A^i^^iVth , (16) 

with A^th from (15) and where Ni^c is defined for the /nl = 1 local model: 



1,2 ,1oc(/nl=1)2 

h 
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Of course, for the local model, the quantities are identical = /^l)- However, when we quote constraints 
on other models we will use Fnl — making self-evident the comparable nature of this quantity — while also 
noting the previously used in the literature. 

The problem with /^^ is that it derives from a somewhat arbitrary normalisation of the primordial bispec- 
trum B^^{ki,k2,k^) which bears little relation to the observable CMB bispectrum signal. The convention 
has been to assume a nearly scale-invariant shape function S{ki, k2, ks) and then to normalise it such that 
S^^{k,k,k) = 1, that is, at a single point; this becomes /nl = 1 case for the model under study. This 
definition ignores the behaviour away from the equilateral value k=k\=k2=k^. For example, S rises from 
a central minimum in the local model and falls from a maximum in the equilateral model; hence, the huge 
disparities between their /nl constraints, e.g. A/^^*^'' « 7A/^£. This definition also does not apply to 
non-scaling models. The alternative to base the non-Gaussianity measure on the actually observable CMB 
bispectrum bf^ij^, as above in (16), docs accommodates non-scale invariant models, such as feature models. 
It also covers bispectra induced by late-time processes like gravitational lensing and cosmic strings. For 
models which are not scale-invariant it should be quoted with the observational cut-off /max- The normalisa- 
tion for a particular model can be easily forecast using the primordial B^^{ki, k2, k^) without the need 
for accurate CMB calculations of bf^i^i^ in (8); primordial shape autocorrelators just need to be compared 
with the local shape as demonstrated in ref. [5] . 
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III. SEPARABLE MODE EXPANSIONS 



When analysing the CMB bispectrum bij^i^, we are restricted to a tetrahedral domain of multipole 
triples {hhh} satisfying both a triangle condition and a limit given by the maximum resolution Imax of 
the experiment. This three-dimensional domain Vj- of allowed multipoles is illustrated in fig. 2 and it is 
explicitly defined by 

Resolution: h,l2,h < ^max , h,l2,h , 

Triangle condition: h < h + h for h > h, h, + cyclic perms. , (18) 
Parity condition: li + I2 + I3 = 2n , n G N . 

The multipole domain is denoted a 'tetrapyd' because it arises from the union of a regular tetrahedron 
from the origin out to the plane h + I2 + h ^ 2/niax and a triangular pyramid constructed from the corner 
of the cube taking in the remaining multipole values out to k < Zmax- Summed bispectrum expressions 
such as (15) indicate that we must define a weight function wij^i^ on the tetrapyd domain in terms of the 
geometrical factor hi^i^is, that is, 

Whhh = hli^h ■ (19) 

This is a nearly constant function on cross sections defined by h + I2 + h = const, except very near the 
tetrahedral boundaries where it is still bounded, and a useful and accurate continuum limit w{li,l2,h) is 
given in [1]. In order to eliminate an /~^/^ scaling in the bispectrum estimator functions, we usually exploit 
the freedom to divide by a separable function and to employ instead the weight 

Ws{li,l2,h) = where vi = {21 + l)^/*" . (20) 

Vj Vj Vi 

tl t2 ^3 

We can then define an inner product of two functions f{lij2jh)j dih^h^h) on the tetrapyd domain (18) 
through 

if, 9) = ^sih,l2,k)fih,l2,h)9{h,l2,h)- (21) 

Given that calculations generally deal with smooth functions /, g, w, v, we can use a variety of schemes to 
speed up this summation (effectively an integration). 

Our goal is to represent the observed CMB bispectrum estimator functions, such as those in (12) and 
(15), on the multipole domain (18) using a separable mode expansion. 



hM,=Y.<Qn{hMM). (22) 



where the Qn are basis functions constructed from symmetrised polynomial products 

Qn{h, h, k) = liQpih) Qrih) Qsik) + Qrih) Qpih) Qsik) + cyclic perms in prs] 

= QipQrqs} with n -H- {prs} , (23) 

with the qp{l) defined below. Here, the six permutations of the polynomial products which we denote 
as {prs} reflect the underlying symmetries of the bispectrum bij^i^ . For convenience, we deflne a one- 
to-one mapping n -H- {prs} ordering the permuted triple indices into a single list labelled by n G N. 
Alternative 'slicing' and 'distance' orderings were presented in ref. [1], but the results presented here are 
robust to this change. However, we shall quote explicit coefficients resulting from distance ordering 

(i.e. n{li, 12,13) < '^'(^1) ^2' ^3) implies Zf + + /| < I'l^ + 1'2 + ^'3 ^'^'^ iii the instance of two modes being 
equidistant the one with most equal li takes precedence). 
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Figure 2: Observational domain (18) for the CMB bispectrum bi-^i^i^. Allowed multipole values {h, I2, I3) lie inside the shaded 
'tetrapyd' region, satisfying both the triangle condition and I < L = /max- 

We choose to define the tetrahedral qp{l) polynomials analogously to Legendre polynomials P„ by requiring 
them to be self-orthogonal with respect to the inner product (21), 

{qpih), qrih)) = Spr , (24) 

with the first few polynomials given by go = 0.074, qi = 0.30(-0.61 + 0, Q2 = 1.2(0.26 - 1.1 Z -h^^) etc. More 
precise expressions and generating functions are given in ref. [1]. As products, the qp only confer partial 
orthogonality on the 3D basis functions Q„, but their use is vital for other reasons, given their bounded 
and near scale- invariant behaviour. 
While the product basis functions Q„ are independent and separable, they are not orthogonal in general 



{Qn, Qp) = Inp / ^np ) 



(25) 



so it is very useful to construct a related set of orthonormal mode functions TZn using Gram-Schmidt 
orthogonalisation such that 

{TZn, Up) = 5np ■ (26) 
Working up to a given order N, the two sets of mode functions are related through 



TZn = '^XmpQp for n,p<N, 

p=0 



(27) 



where Xmp is a lower triangular matrix with 



N 



(A 



Jnp 



{Qn^Tlp) and (7 ^)np = ^Yl{X^)nrKp ■ 



(28) 



Knowing A„p allows us to easily systematically evaluate the expansion coefficients in (22) directly from the 
inner product 



N 



yielding = J]](A^ 

p=0 



(29) 
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Indeed, it is more convenient to present our final bispectrum results in the orthonormal TZn basis, 

blMs=J2^nT^n (30) 

because their orthonormality (26) implies a version of Parseval's theorem. Here, wc note that the expansion 
(30) presumes a spectrum normalised as in (16) to have Fnl = 1, that is, with A'" such that Ylm^n^ ~ 
in the estimator (12). 

To summarise, the Qn{h-,h, hY^ are independent separable basis functions built out of the permutations 
of simple products of the polynomials qp{l), which are well-behaved and bounded over the tetrapyd. The 
Qn's in their easily separable form are employed directly in the bispectrum estimator. However, it is more 
straightforward to present results and to use the inner product (21) with the transformed Tin eigenmodes 
because they are orthonormal; a simple matrix expression (29) relates the expansion coefficients a® and 
using the two sets of basis functions. 



IV. RECONSTRUCTING THE CMB BISPECTRUM 



Now consider the implications of substituting the mode expansion (22) into the estimator (12), while 
exploiting the separability of the Gaunt integral (5) , 



li,mi U't^prs 



^2m2(n)^imi(n) Yi^rnAn) 



I /i,mi 




^ [ [%(n)M,(n)%(n)-6(Mg(n)M,G(n))i^^^ 



(31) 
(32) 



Here, the Mp(n) represent versions of the original CMB map filtered with the polynomial Qp with the 
separated weight function {vl^/C^)'~^, that is. 



(33) 



The maps M^(n) incorporate the same mask and a realistic model of the inhomogeneous instrument noise; 

a large ensemble of these maps, calculated from Gaussian simulations, are used in the averaged linear term 
in the estimator (31), allowing for the subtraction of these important effects. Defining the integral over 
these convolved product maps as cubic and linear terms respectively 



= I d2n%(n)M,(n)M,}(n), 



(34) 



the estimator (12) reduces to a simple sum over the mode coefficients 



(35) 
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where /3« = /3«'="'^ - pf'"". 

The estimator sum (35) is straightforward to evaluate, provided the theoretical model coefficients are 
known. It has been separated into a product of three sums over the observational maps (31), followed by a 
single integral over all directions (34). The actual operations entailed in the estimator sum are only O(Z^), 
so these late-time methods are extremely rapid for direct data analysis and for obtaining variances from 
map simulations. However, we note that the preparatory 'one-off' calculations setting up the orthonormal 
eigenmodes and theoretical CMB bispectra are of order 0{l^). We emphasise that the utility of this approach 
depends on a fairly rapidly convergent expansion for the theoretical bispectrum under study (as indicated 
for almost all models studied to date [5]) and the fact that we have constructed a complete set of orthonormal 
eigenmodes on the observed multipole domain (18). 

There is potentially much more information in the observed (3^ coefficients than just the estimator sum 
(35) which only yields /nl for a given theoretical model. Following the steps above in (31), it is easy to show 
(see Appendix) that the expectation value for for an ensemble of maps with a given CMB bispectrum 
(expanded in modes a'^ with amplitude .Fnl) is 

(/3«) = 5^ FNLaS(Q„, Qp) = Fnl Y1 ^nlnp , (36) 

p p 
so that the averaged estimator (35) becomes 

(S) = ^Fnl E E < ^"f = i^^NL 5^ af = Fnl , (37) 

n p n 

where wc have used (29) in transforming to the TZn basis. (Here wc note again that in this basis N"^ = 
Yln^n"^-) Equivalently, then, in the orthonormal frame we have the simple result 

(^^)=i^NL<, (38) 

that is, we expect the best fit 13^ coefficients for a particular realization to be the a^'s themselves (given 
a sufficiently large signal). Assuming that we can extract the coefficients with sufficient significance 
from a particular experiment, this means that we can directly reconstruct the CMB bispectrum using the 
expansion (30). 



V. VALIDATION USING SIMULATED MAPS 

In our previous paper [1], we have validated this mode expansion approach with simulated maps in a 
WMAP realistic context, accurately obtaining the input /nl as well as a fairly good recovery of the actual 
bispectrum coefficients aj^. In that case, we chose to use the equilateral model because the noise analysis 
requirements were not as stringent. Here, however, we wish to undertake a comprehensive reconstruction 
of the bispectrum, so we have ensured an accurate and robust implementation of the full noise and mask 
subtraction, incorporated in the linear correction term to the bispectrum estimator (12). 

As well as the equilateral model, the other canonical test model we have used to validate our methods is 
the local model, characterised by the shape function (7): 

The local model (39) is dominated by signal from squeezed triangles, e.g. ki <C fa, fa. This behaviour 
is strongly reflected in the CMB bispectrum bij^i^ where the impact of the transfer functions A;(/c) is to 
impose a series of acoustic peaks on the underlying local shape. This dominant signal along the edges 
of the tetrapyd for the local model provides a rigorous test for subtracting the inhomogeneous noise and 
masking effects. This is because these effects also exhibit an overall local shape (as we will discuss in a 
future publication [15]). 
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Figure 3: Recovered spectral coefficients (22) from a single map simulation for a local model with /nl ~ 100. The original 
decomposition coefficients for the local model are shown for comparison (blue). The coefficients were recovered in a 
WMAP-realistic context using the KQ75 mask with inhomogeneous noise added. Error bars (la) are also shown for each mode 
as estimated from 1000 Gaussian maps; note that this variance is roughly constant for all modes. 




Figure 4: Recovered 3D bispectrum using the mode decomposition method (30) from the /nl = 100 local model map simulation 
used in fig. 3. The left panel represents the original theoretical bispectrum used to construct map realisations. The right panel 
represents the recovery in a WMAP-realistic context using the spectrum shown in fig. 3. The main features of the bispectrum 
are evident, including the primary acoustic peak and the high signal in the squeezed states along the tetrahedron edges. 

The a^^ coefficients for the local bispectrum bf^f^i^ are illustrated in fig. 3, having been found using the 
robust CMB bispectrum calculation methods described in ref. [5] (which achieve better than 99% correlation 
with the exact results). With only rimax = 31 Q„ modes it was possible to achieve a 94% correlation between 
the partial sum (22) and the original bi^i^i^ for all the models studied (usually greater than 98%). These 
coefficients were then used to create ensembles of local maps with /nl = 100 (resolution /max = 500) using 
the separable mode map simulation methods described in ref. [1]. As a further test independent of modal 
map-making, the local map simulations from ref. [16] were also used to verify the main conclusions. In 
addition, inhomogeneous noise obtained by coadding WMAP V and W channels was included together with 
the use of a KQ75 sky mask, just as in the original non-Gaussian analysis of WMAP5 data. 

The efficacy of the modal estimator (35) in recovering the correct /nl = 100 from the local maps is 
illustrated in fig. 3, with the ensemble average /nl = 93 it 32.5 found to be in good agreement. We expect 
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the value to be slightly low as we are 94% correlated (the maps have been created using an exact method 
rather than from the modes to ensure robustness). This confirmed the results for the equilateral model 
studied in ref. [1], where it was also shown that Gaussian maps give unbiased results around /nl = 0. The 
recovery of the local bispectrum mode coefficients (34) also proved to be remarkably efficient as illustrated 
in fig. 3 for a typical spectrum obtained from a single map realization. The dominant local modes are 
clearly identifiable above the noise (for a signal of this 3a significance), with the typical variance obtained 
from Gaussian maps also shown. The three-dimensional reconstruction for the local model bispectrum is 
illustrated on the tetrapyd domain in fig. 4 and is comparable with the original bispectrum. The dominant 
features are recovered, including the primary acoustic peak at li ^ I2 ^ I3 ^ 200 and the strong signal for 
the squeezed triangles along the edges of the tetrahedron where one oi k ~ 0. Comparing with results for 
the equilateral model in ref. [1], it is clear that for a measurement of 3a significance or more, we should be 
able to distinguish between families of models which are weakly correlated, such as local and equilateral. 



VI. THE WMAP BISPECTRUM 



We now move on to apply the mode decomposition techniques described and validated in the previous 
sections to the analysis of WMAP5 data. Our aim, first, will be to estimate /nl arising from different 
primordial shapes, some as yet unconstrained in the literature (such as the feature models of section VIII and 
the flattened models of section VII D). Secondly, we aim to provide a full reconstruction of the bispectrum 
from the data, using the same pipeline shown to recover local and equilateral bispectra from simulated 
data. The main emphasis of this work is obtaining fast and accurate convergence for many different shapes, 
rather than a fully optimised estimation. The analysis presented here is intended as a proof-of-concept 
for late time modal estimators of non-Gaussianity, gleaning valuable new information from WMAP rather 
achieving a maximal extraction. For this reason our study has a number of limitations, which we enumerate 
here. We do not implement full inverse covariance weighting in the estimator as in (11) [12], but we adopt 
the pseudo-optimal weighting scheme used by the WMAP team for the WMAP 5-year analysis [13]; we use 
multipoles up to iraax = 500, rather than 1000, since the pseudo-optimal /nl error bars tend to saturate 
above that threshold; finally, we work with WMAP 5-year instead of WMAP 7-year data. The reason for 
not using the latest available dataset is not only that this work started well before the 7- year WMAP data 
release, but is also due to the fact that WMAP 5-year data was originally studied with a pseudo-optimal 
weighting approach, thus making a comparison between our results more straightforward. The present work 
represents the initial implementation of this general approach to analysing non-Gaussianity, rather than its 
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Figure 5: Recovered mode coefficients (22) from the WMAP5 coadded V and W maps. Error bars (la) are also shown for 
each mode as estimated from 1000 Gaussian map simulations in WMAP-realistic context. 
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Figure 6: Recovered 3D bispectrum from WMAP5 data showing the result using the reconstructed mode coefficients shown 
in fig. 5 with the partial sum (30). Several density contours (light blue positive and magenta negative) are shown out to 
k < 500. 

completion. 

After coadding the V and W band data (with the same weights as in the WMAP5 analysis), our first 
step was to extract the mode coefficients from the data, following the procedure summarized eqn (33) 
and (34). In our analysis we chose to compute the first n = 31 modes in (35) because this proved sufficient 
to describe almost all theoretical CMB bispectra on the observational domain /max = 500. The resulting 
estimates will be shown in the following sections. As pointed out in (38), by rotating our recovered 
into the orthonormal frame we obtain the best-fit estimate of the actual bispectrum coefficients a^. The 
mode coefficients obtained from the WMAP5 data in this orthonormal frame are plotted in fig. 5. The 
variance is estimated from 1000 Gaussian map simulations, using the pipeline repetitively in the same 
WMAP-realistic context. 

The mode coefficient extraction from the WMAP5 data was straightforward with both the cubic and linear 
terms contributing significantly to the final result. The late-time estimator (12) is sensitive to all forms 
of non-Gaussianity, in contrast to the two or three separable (and oscillating) modes previously extracted 
from the data using primordial estimators. Despite this increased sensitivity, in principle, making the 
method more susceptible to foreground contamination, our results do not appear to have been significantly 
affected after subtraction by the linear term. This has been investigated through extensive testing, including 
increasing mask size, and we will discuss these issues at much greater length in a companion paper [15], 
characterising the mask, noise and other contributions. It is interesting to note here, however, that the 
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Figure 7: Recovered 3D bispectrum from WMAP5 data showing slices through the data at I = h + I2 + I3 — const.. Shoes 
shown are / = 250, 500, 750, 1000, using the same colour scale as fig. 6. 

mode decompositions can also used to characterise anisotropic contributions, such as the inhomogeneous 
noise (and possibly other contaminants). We will show quantitatively how the action of the linear term 
essentially projects out these spurious bispectrum directions from the cubic term, with the local shape being 
the most affected (as noted originally in ref. [11]). 

The extracted mode coefficients (3^ from fig. 5 can be used to reconstruct the full 3D WMAP bispectrum 
using (30). The result of this partial sum is shown in fig. 6, together with a series of transverse slices 
through the bispectrum shown in fig. 7. Visually the WMAP bispectrum bears a qualitative similarity to 
the local CMB bispectrum already used for pipeline validation, illustrated in fig. 4. As we shall discuss, there 
appears to be some local signal emerging from the WMAP data, but the periodicity of the other features 
does not match well with scale-invariant primordial models. Nevertheless, the orthonormal mode coefficients 
plotted in fig. 5 do not individually show significant deviations away from Gaussianity, presuming the 
accuracy of the nearly constant mode variances which are also plotted. We note at the outset, therefore, 
that the WMAP bispectrum shown in fig. 6 is likely to be the result of cosmic variance (perhaps with 
some residual local signal left-over from the noise/mask subtraction and or other contamination). As well 
as constraining specific theoretical models, we shall test the assumption of Gaussianity more generally in 
section IX by considering a measure of the total integrated bispectrum obtained from the squared coefficients 
(3:^^. In the near future, using the full WMAP7 data set and smaller variances we will expand the scope of 
our mode exploration, including principal component analysis and other statistical approaches [15]. 

Before obtaining specific new constraints, we emphasise again that the extraction of the mode coefficients 
provides a completely model-independent assessment of the three-point correlation function. The ap- 
proach provides far more information than that contained in a simple /nl amplitude parameter extraction 
for particular models. Although obvious deviations from Gaussianity are not apparent from this limited 
WMAP5 analysis (i.e. pseudo-optimal error bars and /max = 500), there remains considerable potential with 
new data sets. For Planck, the sensitivity to primordial non-Gaussianity will improve by up to an order of 
magnitude and so the error bars in fig. 5 will shrink dramatically. The prospects for detection of a large 
NG signal remain completely open. 

VII. CONSTRAINTS ON NEARLY SCALE-INVARIANT MODELS 

Constraints on the bispectrum to date have been for scale-invariant models of separable form, primarily 
on the local and equilateral models, discussed previously. There has been significant evolution over time 
for these constraints as both the CMB data and the estimation methodology have improved. However, 
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as table I illustrates (taken from ref. [3]), there is no compelling and confirmed evidence for a significant 
non-Gaussian signal at this stage. Our purpose in this section is to apply our more general mode expansion 
estimator (35) with our WMAP analysis to obtain constraints on a much wider set of scale-invariant models. 
This method can be applied to any model for which there is good convergence with the given nmax modes. 





Local Equilateral 


Pure cubic 


-58 < /nl < 134 Wl, Komatsu et al 2003 [17] 
-54 < /nl < 114 W3, Spergel et al 2007 [18] 


-366 < /nl < 238 Wl, Creminelli et al 2006 [11] 
-256 < /nl < 332 W3, Creminelli et al 2006 [11] 


Pseudo-optimal 


-27 < /nl < 121 Wl, Creminelli et al, 2006 [11] 
-36 < /nl < 100 W3, Creminelli et al 2006 [11] 
27 < /nl < 147 W3, Yadav Wandelt 2008 [14] 
9 < /nl < 129 W3, Smith et al 2009 [12] 
-9 < /nl < HI W5, Komatsu et al 2009 [13] 


-151 < /nl < 253 W5, Komatsu et al 2009 [13] 


Optimal 


12 < /nl < 104 W3, Smith et al 2009 [12] 
-4 < /nl < 80 W5, Smith et al 2009 [12] 
-10 < /nl < 74 W7, Komatsu et al 2010 [2] 


-125 < /nl < 435 W5, Smith et al 2009 [12] 
-254 < /nl < 306 W7, Komatsu et al 2010 [2] 



Table I: Constraints on /nl°',/n'l"' ' obtained by different groups on the one-year (Wl), three-year (W3), five-year (W5), and 
seven-year (W7) WMAP data releases. The estimators employed are the pseudo-optimal (12), the cubic (the same without 
the linear noise term), and the optimal with fuU-covariance weighting (11). All results were in the context of a primordial 
estimator using separable functions to describe the specific model, unlike the general late-time estimator employed here. For 
further details about the estimator methods employed and the significant evolution of these results over time, please refer to 
the review [3]. 




Figure 8: Predicted 3D bispectrum for the constant model up to h < 500. The same thresholds are employed as those shown 
in the WMAP reconstructions in fig. 4 (after an overall rescaling). 
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A. Constant model 

The constant model S{ki,k2,ks) = 1 is the simplest possible primordial shape with triangles of every 
configuration contributing equally, resulting in a CMB bispectrum bijj^ with features due entirely to the 
transfer functions (as we observed for the acoustic peaks shown in fig. 1). The constant model was motivated 
initially by its simplicity with the large-angle analytic solution (9) for the CMB bispectrum [5]. However, 
the constant shape does have other more explicit physical motivation, such as generation during a slowly 
turning trajectory in multifield inflation, denoted quasi-single field inflation [19]. For nearly scale- invariant 
models, the central values for the bispectrum, bm, all have roughly the same profile but with different 
normalisations. The oscillatory properties of the transfer functions create acoustic peaks located at triple 
combinations involving the following multipole values, / ^ 200,500,800,.... To observe the key differences 
between scale invariant models we must study the bispectrum in the plane orthogonal to the {1,1, Z)-direction, 
that is, the directions reflecting changes in the primordial shape functions. 

For the multipole range Zmax < 500 relevant to the present analysis, we have plotted the 3D bispectrum 
in fig. 8. Here, the dominant feature is the primary acoustic peak stretched along the diagonal of the 
tetrapyd, peaking at I = h = I2 = h = 220 and elongated like an extended balloon from I ~ 100 to / ~ 450. 
Evidence for this primary peak would indicate the presence of a primordial and scale-invariant non-Gaussian 
signal, as emphasised in ref. [5] and investigated quantitatively for the local model in ref. [7]. Observing 
the reconstructed WMAP bispectrum shown in fig. 6 there is a central fluctuation at I ~ 140 but it does 
not extend to larger / as would be expected; see the h + I2 + h = 750 slice in fig. 7 (right) corresponding 
to / 250 where the (apparent) WMAP peak has disappeared. If this measured 3D WMAP bispectrum 
were considered to have any statistical significance then it would mitigate against a scale-invariant model, 
motivating the discussion in section VIII. 

A comparison of the mode coefficients a^^°'^^^^ from the constant model CMB bispectrum shown in fig. 9 
indicates little obvious correlation with the WMAP coefficients ^^"^™3,p (^^igQ plotted). Note that the 
constant model mode coefficients are large for the constant offset n = and for n = 3,4,5 reflecting the 
periodicity of the acoustic peak structure (for /^ax = 500), that is, corresponding to the qpqrQs polynomial 
products with prs = {000}, {002}, {111}, {012} (also with related harmonics at lower amplitude with n = 
9, 10, 11). The mode decomposition estimator (35) yields the quantitative constraint 

^nl'* = ]^Y^ ^^wmap = 35 J i 27.4 , (/§r* = 149.4 ± 116.8) , (40) 

n 

where Fnl is the bispectrum parameter normalised relative to the local model (16), while the lower 
case /nl constraint employs the more model-dependent normalisation using the primordial shape func- 
tion S{k,k,k) = 1. The variance here was determined from using 2000 Gaussian simulations in the same 
WMAP-realistic context. It is clear from this result that there is no evidence — given the present precision — 
for a significant constant primordial non-Gaussian signal. 

B. Local model 

The canonical local shape, which we have already introduced in (39), covers a wide range of models 
where the non-Gaussianity is produced by local interactions. These models have their peak signal in 
"squeezed" states where one ki is much smaller than the other two due to non-Gaussianity typically being 
produced on superhorizon scales. Single-field slow-roll inflation is dominated by the local shape, though 
f^l is tiny [20, 21]. The production of large non-Gaussianity during multiple field infiation [22-24] shows 
much greater promise of producing an observable signal through conversion of isocurvature into adiabatic 
perturbations. Large can also be produced in curvaton models [25-27], at the end of inflation from 
reheating mechanisms [28] and also in more exotic scenarios such as (non-local) p-adic inflation [29] and the 
ekpyrotic scenario [30]. For more comprehensive references and recent examples please refer to the review, 
ref. [31]. 
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Figure 9: Comparison between constant model and recovered mode coefficients for the WMAP5 data. Note tfiat the constant 
model incorporates features entirely due to the transfer functions (the acoustic peaks seen in modes n — 3,4,5), which are 
indicators of its primordial origin. 



The distinct mode decomposition of the local model has already been illustrated in fig. 3, together with the 
3D CMB bispectrum in fig. 4 because they were used to validate this estimator methodology in section V. 
There are a number of differences from the constant model expansion reflecting the dominant signal along 
the edges of the tetrahedron, thus favouring higher order polynomials able to describe this localised signal. 
That is, as well as the periodic acoustic peak signal seen in the constant model (n = 3,4,5), the spectrum 
is otherwise dominated by pure modes n = 9, 15, 26 with prs = {003}, {004}, {005}. The expansion is not 
as rapidly convergent but the eigenmode partial sum achieves a 94% correlation by n = 31. 



0.5 




5 10 15 20 25 



Figure 10; Comparison between local model expansion coefficients and recovered modes for the WMAP5 data. Note the 
relatively slow convergence of the local model and the apparent visual correlation of modes. 

To aid comparison with the recovered WMAP bispectrum, we illustrate both in fig. 10. There appears 
to be some correlation between the two sets of data points which is reflected in the result from the mode 
estimator 

F^^l = 54.4 ± 29.4 = 54.4 ± 29.4) . (41) 

This nearly 2a result is consistent with that obtained by other groups in Table I. In particular, it can be 
compared to the raw result of = 59 it 21 obtained by the WMAP7 analysis, before marginalising over 
foregrounds. The similarities between the recovered WMAP bispectrum and the local bispectrum can be 
observed by comparing the 3D bispectrum in fig. 6 and fig. 4 respectively. There are obvious similarities 
around the edges of the tetrahedron where much of the local signal resides. However, we note an important 
precautionary point. Our analysis of the effects of the noise and mask before subtraction from simulations 
in a WMAP-realistic context indicates that these also have a very nearly 'local' shape. We believe the same 
is also true for likely foreground contaminants which we believe also contribute to the cubic term in the 
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estimator (12); this is indicated also be the significant effect of marginaUsation over foregrounds on /nl 
in the WMAP7 analysis. Our late-time analysis is obviously susceptible to all sources of non-Gaussianity, 
unlike the filtered primordial estimator searching for just a couple of modes. While such contaminants 
appear to have been largely removed from our analysis by the linear term (given their anisotropic nature 
and local shape), the local constraint (41) is clearly more susceptible to systematic effects. Further detailed 
investigations of these effects and the shape characterisation of noise, mask and contaminants using our 
mode expansions is the subject of a follow-up publication [15]. 



C. Equilateral models 



Bispectra dominated by contributions from nearly equilateral triangle configurations, ki ^ k2 ^ are 
produced through the amplification of nonlinear effects around the time modes exit the horizon, which can 
be achieved by modifying kinetic terms, as in the DBI model [32] , or by explicitly adding higher derivative 
terms, such as in K-inflation [see, for example, 33]. For DBI infiation, this leads to non-Gaussianity being 
produced with a shape function of the form [32, 34] 



S{ki,k2,k3) 



kik2k-i{ki +k2 + k^y 



^kf + ^{2kfk,-3kfk])+ {kfk,ki-ikfk]ki) 



■ (42) 



This shape is illustrated in fig. 11, together with ghost inflation [35] and a third distinct single field equilateral 
shape found in a general analysis of such models [33]. Note that the generic equilateral shapes are not 
separable, but have been approximated to date using a separable ansatz commonly called the 'equilateral 
model' [11]: 



S^^^\ki,k2,k3) 



1 {k2 + k3- ki){k3 + fci - k2){ki + k2- fcs) 
N kik2k3 



(43) 



Despite the apparent visual differences between these primordial shapes, particularly near the edges of the 
tetrahedral domain, the resulting CMB bispectra share at least a 95% or greater correlation ([see 5]). The 
CMB mode decomposition for these models is illustrated in fig. 12, showing very similar behaviour to the 
constant model also dominated by the acoustic peak coefficients n = 3, 4, 5. The resulting constraints from 
the modal estimator are: 



Equilateral: 
DBI: 
Ghost: 
Single: 



Fnl 



25.1 ± 26.4 

26.7 ± 26.5 
22.0 ±26.3 

28.8 ±26.6 



(/nl = 143.5 ± 151.2) , 
(/nl = 146.0 ± 144.5) , 
(/nl = 138.7 ±165.4), 
(/nl = 142.1 ± 131.3) . 



(44) 
(45) 
(46) 
(47) 




Figure 11: The sliape function of models in the equilateral class which from left to right are DBI inflation, ghost inflation and 
the remaining single field inflation model. 
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Here, the local i^NL normalisation (16) yields much more consistent variances between models within the 
equilateral family than /nl (as well as values comparable to local and other models). Note that there are 
up to 30% variations between the central values of these Fnl constraints despite the strong correlations 
between these bispectra; this is attributable to the different behaviour near the edges and faces where 
much of the apparent WMAP signal is localised. These results are consistent with the evolving constraints 
obtained in the literature to date, as shown in Table I. 

Finally, we consider a separable 'orthogonal' shape 5'°'^*'^°^ which is a constructed from a linear combination 
of the constant and equilateral shape functions 5™**^°^ oc 5"^^"'^ — 2/3 (see [12, 36]). The constraint from 
the mode estimator (35) then becomes 

^ortho ^ 3 ^ 27.3 , {f^t^° = -79.4 ± 133.3) , (48) 
which is a less negative result than the latest WMAP7 limit f^i^° = —199 it 104, but it remains consistent. 



D. Flattened model 



It is possible to consider inflationary vacuum states which are more general than the Bunch-Davies 
vacuum, such as an excited Gaussian (and Hadamard) state [37, see also discussions in Chen et al. 33, 
Meerburg et al. 36] . Observations of non-Gaussianity in this case might provide insight into trans-Planckian 
physics. The proposed non-separable shape for the bispectrum is 

The bispectrum contribution from early times is dominated by flattened triangles, with e.g. ~ ki + k2, 
and for a small sound speed <C 1 can be large. Unfortunately, as the divergent analytic approximation 




-Constant 

"Single 

-DBI 

"Equilateral 
"Ghost 




Figure 12: Equilateral model expansion coefficients compared between models (top panel) and compared with WMAP5 
results (lower panel). 
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breaks down at the boundary of the allowed tetrahedron, some form of cut-off must be imposed, as shown 
for the smoothed shape in fig. 13 where an edge truncation has been imposed together with a mild Gaussian 
filter. This leads to a degree of predictive uncertainty, but the regularisation scheme ensures the primary 
signal is well-localised on the tetrahedral faces and is quite distinct from other separable shapes investigated 
to date (refer to ref. [5] for the specific details). 

The resulting CMB spectrum reflects this behaviour with the dominant signal residing near the tetrahedral 
faces as shown in fig. 13. Figure 14 shows the flat model mode coefficients, which like the local model are 
only slowly convergent. Comparing the flat model with the coefficients obtained from WMAP, the mode 
estimator yields the new constraint: 

Fnl = 35.4 ± 29.2 (/nl = 18.1 ± 14.9) . (50) 

Despite the apparent visual similarity between the flat and WMAP CMB bispectra (figs 13 and 6) the 
present analysis does not reveal a significant correlation. 

E. Warm model 

Finally, we consider warm inflation scenarios, that is, nearly scale-invariant models in which dissipative 
effects play a dynamical role, because these also may produce significant non-Gaussianity [38] (for a review 
see [39]). Contributions are again dominated by squeezed configurations but with a different more complex 
shape possessing a sign flip as the corner is approached, essentially making the signal orthogonal to the 
local shape. It can be shown that this makes the warm and local shapes essentially orthogonal with only a 
33% correlation (see ref. [5] where the shape function and CMB bispectra are discussed). As with the flat 
model, uncertainties remain as to the validity of the approximations made as the corners and edges of the 
tetrapyd are approached. Comparison of the predicted warm bispectrum coefficients (3^ with the WMAP 
data through the modal estimator (35) yields the constraint 

F]^^™ = 10.3 ± 27.2 (/i^f = 47.4 ± 125.4) . (51) 

A previous WMAP3 warm inflation analysis obtained a lower central value /nl'^™ = — 169± 103 [40] which 
is marginally consistent with (51) at the 95% confidence level. Probably the most significant difference 
is that the previous analysis did not include a linear term in the estimator (12) to account for noise and 
masking effects; these corrections are significant here as for the edge-dominated local model. 




Figure 13: Flattened model: smoothed primordial shape function (left) and three-dimensional CMB bispectrum (right) for the 
flattened model. 
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Figure 14: Flat model mode coefficients compared to WMAP5 mode coefficients. 



VIII. IMPLICATIONS FOR NON-SCALING FEATURE MODELS 



It is possible to produce non-Gaussian signals which are not scale-invariant, such as models with a distinct 
feature in the inflaton potential. These usually take the form of either a step in the potential (models which 
have a long history, see e.g. ref. [41]) or those with a small oscillation superimposed onto the potential 
(which have become more popular recently, see e.g. ref. [42]. Two analytic forms for the resulting three 
point functions have been presented in ref. Chen et al. [41] with the expression we will analyse here taking 




Figure 15: Feature model coefficients plotted in two-dimensions by mode number n and as function of phase tp with I* — 400 
(top panel) and as a function of scale I* with (p = (lower panel). Note how the characteristic n = 3,4, 5 primordial acoustic 
peak signature is affected (compare with fig. 9. 
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Scale 

Phase^~"~~---,,,_^ 


150 


200 


250 


300 


400 


500 


600 


700 





57 (30) 


-52 (33) 


-25 (32) 


1(30) 


1(27) 


8 (26) 


18 (25) 


23 (25) 


7r/8 


67 (36) 


-26 (27) 


-36 (30) 


-6 (25) 


-4 (26) 


-2 (27) 


12 (26) 


20 (25) 


7r/4 


68 (42) 


-10 (29) 


-43 (30) 


-11 (21) 


-7 (25) 


-10 (27) 


-1 (28) 


13 (27) 


37I-/8 


49 (46) 


7 (34) 


-42 (32) 


-18 (24) 


-9 (25) 


-14 (26) 


-13 (28) 


-2 (28) 


7r/2 


15 (46) 


32 (41) 


-30 (35) 


-32 (34) 


-10 (25) 


-16 (25) 


-18 (27) 


-14 (28) 


57r/8 


-19 (42) 


63 (46) 


-15 (35) 


-38 (43) 


-11 (25) 


-16 (25) 


-20 (26) 


-20 (27) 


37r/4 


-39 (35) 


87 (48) 


(35) 


-25 (41) 


-11 (26) 


-15 (25) 


-21 (25) 


-23 (26) 


77r/8 


-48 (30) 


81 (43) 


13 (34) 


-11 (35) 


-7 (27) 


-13 (25) 


-20 (25) 


-23 (25) 



Table II: Limits for a selection of feature models in the form Fnl (StDev). 



the form 

S^^'^^h, k2, h) = 1 sin (2^^^^^^ + ^) , (52) 

where k* is the associated with the physical scale of the feature in question and $ is an arbitrary phase 
factor. The alternative form with a logarithmic momentum dependence in the sin argument can be shown 
to be closely correlated with the simpler form (52), certainly on the present domain of study /max = 500. 
Previously, we studied the shape and CMB bispectrum for a particular feature model (with k* «i T/tq and 
I* ~ 400), showing that its non-scaling behaviour made it essentially independent of all the other shapes 
[5]. Such models can have starkly contrasting CMB bispectra as illustrated in fig. 18, disrupting the usual 
pattern of acoustic peaks which switch from correlation to anticorrelation on multipole scales 1* . Clearly, 
scale dependent feature models form a distinct category of bispectra beyond the equilateral, local, warm 
and flat families, so searches within WMAP and future data sets are well-motivated. 

For the present WMAP5 analysis, we have studied the primordial feature shape (52) over a wide range of 
for which the CMB bispectra that we obtained could be accurately described by our n = 31 eigenmodes, that 
is, for which we could obtain > 90% convergence to bff\ for the partial sum (22). This restricted the scale 
parameters in (52) to the range /* > 150, so wc studied Values /* = 150, 200, 250, 300, 400, 500, 600, 700. 
For larger values /* > 700 the models became highly correlated with the constant model given that Imax = 
500. No such restriction applied to the phase which was studied for each I* over the full domain < $ < 27r 
in 7r/8 steps (noting that models separated by vr are merely anticorrelated). This entailed considerable 
computational effort calculating 64 distinct CMB bispectra at high accuracy using the robust methods 
previously described elsewhere [43]. The mode coefficients for the I* = 400 model are illustrated for the 
different phases in fig. 15, demonstrating how the characteristic acoustic peak signal in n = 3, 4, 5 can 
be modified (compare the constant model fig. 9). The strong dependence of the mode coefficients on the 
different multipole scales I* (at fixed phase $ = 0) are shown in fig. 15. 

Results from the modal estimator for all the feature models investigated are provided in Table II. Note 
that the constraints are given in terms of the normalised quantity Fnl defined in (16), since there is no 
simple generalisation of the primordial normalisation used for /nl without scale-invariance. As before, 
the variances (given in parentheses) are those obtained for the same set of models from 1000 Gaussian 
simulations. The results are illustrated graphically in fig. 16 showing the relative significance of the central 
Fkl values relative to the standard deviation. The result with the highest significance is that for the feature 
model with 1* = 150 and zero phase which achieves a 1.9a significance. The 3D bispectrum for this model 
is shown in fig. 18 demonstrating how such models can reproduce the apparent scale-dependence observed 
in the WMAP bispectrum (sec fig. 6). However, we note that this model is close to the resolution limit set 
by the eigenmodes deployed (like the other cases of higher significance). The results over the full domain 
of feature models investigated remain consistent with the Gaussian hypothesis with no significant detection 
found on the WMAP domain for I < 500. 



22 




Phase (X 71/ 8) ^ Scale 



Figure 16: Significance of feature model bispectra _Fnl/A^nl using WMAP data with the modal estimator (35). This is plotted 
as a function of the multipole scale l* and the phase of feature models given by (52). 
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Figure 17: Best fit feature model coefficients {1* = 150, </> = 0) compared to WMAP5 mode coefficients. 
IX. TOWARDS A MEASURE OF THE TOTAL INTEGRATED BISPECTRUM Fnl 



Our focus in this paper has been on recovering the observed bispectrum bi-^i^i^ which contains more 
information than constraints for particular models. We can also consider squaring this quantity and 
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summing over all multipoles to obtain a total integrated nonlinearity parameter defined by 

2 



Substituting our mode decomposition (30) into the expression for the integrated bispectrum we can find 
the leading order contribution from the three-point correlator to be^ 



tc2 



- 7^ E E m n,) ^ . (55) 

Unfortunately, however, the expectation value (-F^l) contains much more than the just contributions 
from the three-point correlator. There is necessarily a contribution from products of two-point correlators. 



As an unambiguous signature of a significant bispectrum we sfiould compare -Fj^l witli tlie skewness 71 wiiicii is given by [44] 

71 = ^(^(ft))'^ = ■ (54) 

fn principle, the skewness can conspire to vanish even with a non-zero bispectrum bi-^i^i^ because it is not positive definite, in contrast 
to the bispectrum contribution to -Fnl- 
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reflecting both cosmic variance and instrument noise and there are also higher-order contributions, which 
are derived in the Appendix. The leading order terms are found to be 



1 



NL 



Ti2 



(56) 



where the first term represents the underlying variance, the second is the integrated bispectrum we seek, 
and the third term arises from contributions from the four-point correlator or trispectrum t1^i^{L). 

The non-Gaussian corrections from the bispectrum and trispectrum (and possibly higher correlators) to 
Fj^L easily distinguishable. Nevertheless, we know the Gaussian expectation so it can still be used 

to measure any deviation from Gaussianity, by determining the amplitude and nature of any differences. 
To prove the utility of this statistic we have applied it to both Gaussian and local simulations. As we have 
shown in a previous paper [44] our local map simulations contain a negligible trispectrum and so we can 
use the approximation 

^NL ^ <^Nl'> + T^i^NL E ' (57) 



loc 



a. 



tc2 



where we have defined {F^^^ to be the value recovered from Gaussian simulations. As iV^ = ^ 
and we normalise all models to have N = N^°^^ we can then attempt to recover without assuming any 
particular from for 6i^. The recovered Fnl will then be given by the following formula 



NL 



F2 



(58) 



We estimated (^nl^) fro^i 1000 Gaussian simulations and then calculated from 100 simulations with 
various input Fnl • The results are presented in table III and are plotted as cumulative sum of squared mode 
coefficients in figure 19. The encouraging results show that the statistic is recovering the input Jml to a 
reasonable degree of accuracy, although with a slight tendency to under-estimation. The biggest surprise is 
that the error bars for this general bispectrum estimator only increase by a factor of 50% from those when 
the particular local form is assumed explicitly. We note, however, that (57) gives rise to a x^-distribution, 
so we have to take care in assuming Gaussianity for small rzmax- Again, we will further explore the utility 
of such general modal statistics elsewhere [15]. 



Input /jvL 


Mean 


StDev 


Null 


50 


55.86 


48.53 


36 


75 


65.39 


46.95 


20 


100 


95.28 


41.61 


8 


150 


138.8 


44.87 


3 


200 


187.39 


41.55 






Table III: i^Nif as recovered from 100 simulated local maps. Null refers to the number of maps in which the recovered Jnl is 
less than the Gaussian expectation value. 

With the efficacy of the Fnl statistic established we have also applied it to the WMAP5 data. This 
yields the unexpected result that Fnl obtained from the WMAP5 data is less than that which we we 
would expect from a typical Gaussian map by slightly over la (see the cumulative sum in fig. 20). This 
is somewhat surprising because one would expect a late-time estimator to be susceptible to foregrounds or 
other contamination, but the deviation remains statistically insignificant. In principle, this result could be 
due to a large negative trispectrum or higher order contribution (see Appendix). However, neglecting this 
possibility, the result shown in fig. 20 indicates that there is no significant contribution to the bispectrum 
from the first 31 eigenmodes. This would constrain virtually all smooth scale invariant shapes, as well as 
the feature models we have surveyed. The only remaining possibility for a bispectrum detection (at the 
present precision) would then be for oscillatory models with sufficiently high frequencies or bispectra with 
particularly sharp, or localised, features (i.e. those which require n > 31 for an accurate description). We 
have good evidence, therefore, for the null hypothesis that we live in a Gaussian universe. 



25 



500 
450 
400 
350 
300 
250 




— Gauss 
^f,,=100 

-W=200 ' " " " " 


1 [ 1 

- - - - - 




- - - - /' - - 


- — 


200 
150 




y-^^r^r^rr^.. , 


— .^^^^^ — - 


100 
50 













1 


1 1 



5 10 15 20 25 30 

Mode 

Figure 19: Cumulative sum of mode contributions to the total (55) for the local Jnl ~ 100 (red) and Fnl = 200 (green) 
map simulations compared with Gaussian maps (blue). The la variance is shaded around the mean value obtained from 100 
simulations (1000 simulations for the Gaussian case). 
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Figure 20: Cumulative sum of mode contributions to the total Jnl (55) for the WMAP data compared with Gaussian map 
simulations as in fig. 19. 

X. DISCUSSION AND CONCLUSIONS 

We have implemented and validated separable mode expansions with a general late-time CMB bispectrum 
estimator, using it to investigate a wide range of primordial models with WMAP 5-year data. Notable new 
constraints include those on non-scaling feature models, trans-Planckian (flat) models and warm inflation. 
The results for nearly scale-invariant models are summarised in Table IV, demonstrating consistency with 
previous constraints on equilateral and local models. Note that we adopt a nonlinearity parameter Fnl 
normalised to facilitate direct comparison between the local /nl and any other model. We found no evidence 
for significant deviations from Gaussianity for any specific model (at 95% confidence). Feature models were 
surveyed over a wide range of parameters with periodicities above /* = 150 and over the full domain of 
phase values. Again, no significant bispectrum detection was made, though given the nature of this survey 
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Model 




(./■nl) 


Constant 


35.1 ± 27.4 


(149.4 ± 116.8) 


DBI 


26.7 ±26.5 


(146.0 ± 144.5) 


Equilateral 


25.1 ± 26.4 


(143.5 ± 151.2) 


Flat (Smoothed) 


35.4 ± 29.2 


(18.1 ± 14.9) 


Ghost 


22.0 ± 26.3 


(138.7 ± 165.4) 


Local 


54.4 ± 29.4 


(54.4 ± 29.4) 


Orthogonal 


-16.3 ±27.3 


(-79.4 ± 133.3) 


Single 


28.8 ±26.6 


(142.1 ±131.3) 


Warm 


24.2 ± 27.3 


(94.7 ± 106.8) 



Table IV: Limits for all known scale invariant models 



some models provide a better a posteriori fit to the data than others. 

As we have emphasised throughout, more information can be extracted from the mode decomposition of 
the data than a few Fnl's for specific models. Given that we have constructed a complete orthonormal 
basis TZn we can use the mode coefficients to directly reconstruct the full CMB bispectrum using the 
partial sum (30). We plotted the result for WMAP5 in fig. 6 which, despite its low significance, revealed 
interesting qualitative features similar to the local model (4), but without the periodicity expected from 
acoustic peaks. We discussed a positive-definite measure for the total integrated bispectrum constructed 
from the mode coefficients = X^n /^n ^) which was used to recover /nl from map simulations in a model 
independent manner (though with larger variance). For WMAP5 data the integrated Fnl was found to be 
small and again consistent with a Gaussian hypothesis. 

Despite the absence of any convincing evidence for a statistically significant CMB bispectrum in the 
present analysis, many avenues remain open for further investigation using the present methodology. The 
late-time modal estimator (35) can identify any bispectrum whether generated at early times like inflation 
or sourced since decoupling by cosmic strings, gravitational lensing, or second-order gravitational effects. 
Unlike the primordial estimator, the general mode expansion can also be used to characterise noise and 
foregrounds, which need to be identified and subtracted through the linear term in the estimator (12). 
The efficacy of this removal and other validation checks which may affect a residual local signal will be 
published shortly [15]. Finally, we note again that these methods can be pressed much further with existing 
and future data, especially from Planck. The anticipated Planck variance A/nl ~ 5 will substantially 
improve sensitivity to specific bispectrum shapes, leaving significant discovery potential available in the 
near future. We note also that these separable mode techniques have been adapted for general CMB 
trispectrum estimation, in principle, making tractable the investigation of all planar primordial trispectra 
[44] . Analogous methods can also be applied to modal bispectrum extraction for large-scale structure and in 
other contexts. For the time being, however, this general bispectrum survey uncovers no significant evidence 
of non-Gaussianity which would undermine the standard predictions of the simplest models of inflation. 

XI. ACKNOWLEDGEMENTS 

We arc very grateful for many informative discussions with Donough Regan, Xingang Chen, Anthony 
Challinor and Alessandro Renzi. Simulations were performed on the COSMOS supercomputer (an Altix 
4700) which is funded by STFC, HEFCE and SGI. We are particularly grateful for computer support 
from Audrey KaHazin. JRF, ML and EPS were supported by STFC grant ST/F002998/1 and the Centre 
for Theoretical Cosmology. EPS is grateful for the hospitality of the Arnold Sommerfeld Centre and the 
Universe Excellence Cluster in Munich. 



27 



'3 



XII. APPENDIX - HIGHER-ORDER CONTRIBUTIONS TO 

Consider the expectation value of the quantity {F^i) = J^ni^n"^) which we defined in (53) and discussed 
in section IX. Our definition of the bispectrum mode coefficients in (34) can be expanded using the 
expressions in (31) to the expficit form 

5k hij^i^TZn{ll,l2,h) Sr^ ( h h h\ /c-n\ 
Pn = 2^ ^ ^ 7. ahm^ahm^O-hrm ■ (59) 

^ vi^vi^vi.^yCi.Ci^Ci^ ^ \ mi 1712 m3 J 

It is instructive at this point to repeat the derivation of the expectation value of the coefficient {/3J^) for an 
ensemble of universes with both a given Fnl as in (16) and a given theoretical bispectrum ^;^/^'^^~^^ (for 

further details see ref. [1]). We describe the theoretical bispectrum b^^^^i^^~ by the orthonormal mode 
expansion coefficients as in (30) or, equivalently, the separable modes as in (22), distinguishing it 
from the observed bispectrum hj^i^ described by or . Substituting the separable mode expansion of 
the bispectrum (22) into the expression for the recovered coefficient (34) wc find 

ip-) = /d2n(M{,(n)M,(n)M,}(n)) = ^ (g^^^^ f (6,^,^,^) (60) 

= Fnl"^ WsQnQp = Ff^h'^lnpap . (62) 

P hhh P 

^Prom (29) we can transform this into the orthonormal basis TZn to find the simple result 

Wn) = FNLa^ ■ (63) 

Here we have ignored the linear term in (34) because its expectation value vanishes. 

Now the expectation value of the square of this coefficient necessarily involves the six-point function, 
because the expression takes the form 

imi (i'hm^ ^hms (^Um4 ^(sms ^hme, 

^ \^ mi m2 ms J ^1714 me J 

Here we note that we can include a cut-sky, as well as noise and beam effects, as we did previously in (14) 
and (13) respectively. 

The expectation value of the six-point function has a variety of non-vanishing contributions from combi- 



(61) 
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nations of lower order correlators, which become (after summing over equivalent permutations): 



Ia ^5 ^6 



-1112 ^rri s 



C W ^2 ^3 ^4 ^5 ^6 pi A 



6 



I QH ^1 ^2 ^4 ^3 ^5 le 



m\m2mi^ m^m^me 

+ GSi^lJrm -m2 Cli X!*^"-"-^ 
LM 

i^mi —mi 

LM 

+ higher order terms , 



2t4 ^ishle 
M ^ 

M I 



h h L 
ms m4 M 

h h L 
m2 ms M 



h h L 



1715 rriQ 



L 

-M 




(65) 



where the CMB trispectrum T^^^'^{L) is reviewed and discussed at some length in a recent companion paper 

[44]. 

Consecutively labelling the ith term in the six-point expression (65) above, we now evaluate each specific 
contribution to the expectation value of -F^l, denoting modes term-by-term as {Pn'^)i- Recall that we have 
orthonormal basis eigenmodes Tin (26) which satisfy 



E 



Ws{li,h, h)'Tln(,h, h, h)T^pih, h, h) — {T^n, 'Tip) — ^np ■ 



(66) 



The first term from (65) is a product of two-point correlators Q in the numerator which cancels with the 
weighting in the denominator of (64) to become simply 

(;9f)i = 6(^„,:^„)=6. (67) 

This is the primary Gaussian noise contribution which cumulatively dominates -P^l, as discussed in section 
IX, and which has been confirmed quantitatively in Gaussian simulations. The second term from (65) also 
consists of products of two-point correlators which divide out, but it does not simply further: 



(^r)2 = 9 

h 



iP xP" iP 

li I3 /s 




(68) 



Despite the difficulty evaluating this expression explicity, it appears that the product of Wigner-3j symbols 
will behave asymptotically as so we can expect the summed product of the Tin eigenfunctions in (68) 
to be significantly suppressed relative to the inner product in (67). For ^max 3> 1, therefore, we expect 
^ (/^n^)2- For small non-Gaussianity, it may be necessary to calculate this term explicitly, though 
it is more costly to evaluate than the usual inner product (21). Alternatively, its effect can be determined 
from Gaussian simulations, which already confirm the clear dominance of the first contribution (67). 

The third contribution from (65) is the straightforward product of the bispectra sought in IX, which from 
(30) simply collapses to 



^^1/2^3 



v} v} 



\/Ci^Ci^Ci. 



--Tln{h,h,h)bi^ 



w, 



:{ll,l2, h)Tln F^h E Q^Tlr, 



p2 /^7i2 



(69) 



where here we have distinguished the recovered hi^i^i^ in the above from that normalised with Fnl = 1 which 
defines the coefficients in (30) . The fourth contribution also arises from products of the bispectrum but 
in the form 

2 

1 



E 



2L + 1 



E 



a 



I1I2 



Ws{ll,l2, L)Tlni{ll,l2, L)Tlp{{li, I2, L) 



(70) 
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However, as for the second term (68), it is clear that the additional weighting (2L + 1) ^ will generically 
suppress this contribution, so that {Pn'^)^ <C (/3^^)3 for I 

max 3> 1. Finally, we consider terms involving 
combinations of the two-point and four-point correlators which simplify considerably because of the following 
identities for summed Wigner-3j symbols: 



E 



m2,m3 



h h h \ ( h h ^4 \ _ ^hh^-i 



mi—nii 



mi 

The fifth term from (65) reduces to 



mi 7712 I \ m2 7714 / 2Zi -|- 1 

. [ ^1 ^1 h \ ^ y2Zm 5/30 ■ (71) 
\ TTii —mi / 



which is a monopole term which can be ignored. The sixth term from (65) can be expressed in the form 
= 9 V — V V(2^i + 1)(^^2 + 1)(2^3 + 1)(2^4 + 1) ( L h l2 \ ( L h lA 

^Pn /6 ^Vl^ \T,V,^V,^Vi,Vi,^Cifiifii,Ci, 1^ ^ 1^ ^ 

X 7^„(L, hM) nn{L, h, k) Tl±{L) , (73) 

This term will contribute at a similar order to the third term with (/3^^)3 ~ if there is a significant 
trispectrum, that is, if S'nl'Tnl 3> 1 (see ref. [44]). Assuming that correlators (unconnected parts) beyond 
fourth-order are negligible, then the primary contributions to F^-^ become 

FL = 6n„.ax + J2 [^NL«n' + {Pn\ + -] ■ (74) 
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